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Second order gauge invariant discretizations 
to the Schrodinger and Pauli equations 

Snorre H. Christiansen* Tore G. Halvorsen^ 


Abstract 

We introduce a numerical method, based on finite elements and lattice 
gauge theory, to compute approximate solutions to Schrodinger and Pauli 
equations. The crucial geometric property of the method is discrete gauge 
invariance. The main new achievement is second order convergence. This 
is proved by iirterpreting the method as defined on gauge potential de¬ 
pendent finite element spaces and providing an analysis of such spaces in 
terms of gauge potential dependent norms on simplices of all dimensions. 


1 Introduction 

The Schrodinger equation is the fundamental equation of non-relativistic quan¬ 
tum mechanics. It couples to an electromagnetic field through an associated 
gauge potential (see for instance chapter 4 in m)- Since there is some freedom 
in choosing a gauge potential for a given electromagnetic field, it is important 
that the solution to the Schrodinger equation transforms, when the gauge po¬ 
tential changes, in such a way that observable quantities, such as energy levels 
and probability densities, remain unchanged. 

This paper is concerned with developing a finite element method with a 
similar property, and it is inspired by lattice gauge theory in the sense of M- 
The new method improves upon our previous works [7] and [S] by yielding a 
higher order of convergence. In order to be more precise, on the problem and 
our new results, we need to introduce some notations. 

We let S denote some spatial domain in which we assume to be bounded, 
convex and polyhedral. 

A gauge potential on S is just a vector field. Given one, called A, we consider 
the covariant gradient, defined on complex valued functions on S by: 

Vau = Vm -I- iAu. (1) 

We will be interested in taking scalar products of complex functions, spinors, 
as well as complexified vectors and one-forms. The C-bilinear scalar products 
will be denoted: 

(u, v) ^ u - V. (2) 
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Hermitian scalar products then take the form; 

(u, v) 1 -^ • V, 

where denotes the complex conjugate of u. 

The adjoint of the operator Va with respect to the hermitian scalar 


products on S is denoted We may remark that: 

=—{V ■ u + iA ■ u). (4) 

We then define the covariant Laplacian: 

Aa = -V:^Va. (5) 

We may expand this expression as follows: 

Aau = Au + 2iA ■ Vu + i(V ■ A)u — (6) 

We may also define the covariant differential operators: 

gradA M = Vam = gradu + iAu, 
curlA M = VA X u = cnAu + iAxu, (7) 

divA M = Va ■ w = div u + iA-u, 

in terms of which we have: 

Aam = divA gradA u. (8) 


We are interested in eigenvalue computations of the following form. Find a 
complex function u on S and € R such that: 

— Aam = Eu. (9) 

We use Dirichlet boundary conditions, that is u\gs = 0. In the last section of 
the paper we consider some extensions of this eigenvalue problem: inclusion of 
a scalar potential and taking into account spin. The latter is done through the 
Pauli equation. However for the remainder of this introduction we stick to (P|). 

The variational formulation is to find u € Hj(S') and Fi in R such that for 
all r; e Hi(S'): 

a{u^v) = E(u,v), (10) 

with: 

(m, ~ J (H) 

and: 

a{u,v) = J{VAu)^ ■ {VAv). (12) 

Depending on the situation we will include or not the dependence of a on A, by 
writing a[A\. 

A crucial property of this eigenvalue problem is gauge invariance. Given a 
scalar field a on S', with real values, we may transform A and u as follows: 

A A' = A — grad a, 
u' = exp(iQ;)tt. 
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(13) 

(14) 


Then we have: 


( 15 ) 


V A'u' = exp(iQ;)VAU. 

We therefore have the invariance properties: 

{u\v') = {u,v), (16) 

a[A'](u' ,v') = a[A\(u,v). (17) 

These invariance properties are related to local conservation of electric charge, 
via Noether’s theorems. 

For the eigenvalue problem, they have the consequence that if: 

— ls.AU = Eu, (18) 

then: 

- Aa'u' = Eu', (19) 

Moreover, concerning the associated probability densities, we notice that they 
are the same for u and u'. That is, for all a; £ S' we have: 

|n'(:r)r = Kx)r (20) 

We wish to construct a numerical method with similar invariance properties. 
Consider regular simplicial meshes 7h of mesh-width h. On Th we have the 
space Xh of complex-valued continuous piecewise affine functions. The standard 
Galerkin method is to find u £ Xh and Ti in R such that for all v £ Xh'. 

a[A\(u,v) = E{u,v). (21) 

The standard results on eigenvalue approximation (see [I]), show that eigenvec¬ 
tors converge with order h in norm and with order h? in norm, whereas 
the eigenvalues converge with order /i^. 

This method is not gaugeinvariant. We suppose that the field A used in 
(1^ is a Whitney one form (defined in section (lO) . or equivalently a Nedelec 
edge element vector field. It is then natural to consider gauge transformations 
for A of the form m with a a Whitney zero form, that is a scalar continuous 
piecewise affine function. However then the gauge trasnformation (11411 maps u 
out of the space Xh- 

A solution is to keep gauge transformations (1131) acting on A, but modify 
the gauge transformations CSl) acting on u, so as to stay within Xh'- 

A^ A' = A — grad a, (22) 

u^ u = n?i exp{ia)u, (23) 

where n?i is the nodal interpolator onto Xh- An interpretation is that one 
modifies just the nodal values by the gauge trasnformation, but stay piecewise 
affine. Then the problem is that we no longer have the invariance of the bilinear 
forms nuini). 

In [7] we proposed a modification of the bilinear forms, inspired by lattice 
gauge theory mm. such that the modified bilinear forms were invariant under 
the discrete gaugetransformations (EH USD, yet stayed close to the original ones. 
In the case of smooth gauge potential we obtained the estimate, for u,v G Xh'. 

\a{u,v) -d{u,v)\ ^/i||M||Hi(s)lkllHi(S)- (24) 
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This guaranteed that the eigenvectors converge at a rate h in H^. However it 
seems that the order of convergence in was just h, and that the order of 
convergence for the eigenvalue is also just h. Moreover our error estimates were 
just valid for meshes for which the discrete maximum principle is true. This 
condition is typically enforced by requiring that dihedral angles be weakly acute. 

In [5] we proposed, among other things, a more elaborate method for such 
eigenvalue problems, which did not require such restrictions on the mesh. This 
was achieved by no longer relying on mass-lumping techniques. However the 
basic estimate is still (El, so that the orders of convergence were still just h. 

The purpose of this paper is to present and analyse a method which is gauge 
invariant, but where second order convergence holds. We also include spin in 
the discussion. Importantly, the theoretical underpinnings of the method, which 
might be of a broader interest, are of a rather new type. 

The underlying idea is different from our previous works. Even if our method 
can be interpreted as a variational crime on Xh, the analysis relies heavily on 
interpreting it as a variational crime on another discretization space. This im¬ 
plicit space, denoted is a discrete space still having one degree of freedom 

per vertex, but where the local shape of the functions is determined by solving 
a local PDE related to the global PDE we are addressing. Explicitely we con¬ 
sider functions u such that for all simplexes T in the mesh, of all dimensions, 
= 0, where A is the average of A on T. For instance we may remark 
that = Xh (affine functions on simplices are characterized by the prop¬ 

erty of being harmonic on all subsimplices). When A ^ 0 we cannot compute 
explicitely the solutions to these local PDFs, except on edges. As it turns out 
this will be enough for defining our numerical method. 

The idea that the discretization space should incorporate the behavior of 
the PDE is not new. Multiscale finite element methods are often based on 
this idea. The method of m can also be interpreted this way, and we share 
with this method that we are especially concerned with the behavior of discrete 
functions on edges, whereas the previously mentioned multiscale methods are 
mostly concerned with the behavior on the maximal simplices (tetrahedra in 
three space dimensions) often using standard finite elements on the skeleton. 
We have previously introduced modified shape functions adapted to convection 
diffusion problems, in the framework of finite element systems mm- It has been 
our hope that the type of analysis we present here might extend to convection 
diffusion equations, but so far this har not been realized. For convection diffusion 
equations one is interested in the regime of vanishing viscosity, whereas in the 
present case we are not considering the highfrequency regime, even though this 
could be interesting in some experimental setups. 

As in the framework of finite element systems, we insist on the recursive 
nature of our discrete functions: they are defined not only on tetrahedra (top¬ 
dimensional simplexes in our mesh) but also on all the subsimplicies of all di¬ 
mensions. While we have previously been much concerned with the algebra 
of this recursive structure [1][5] (for mixed finite elements or, more generally, 
differential forms) we have also analysed stable interpolation operators using 
recursively defined norms nni (e.g. Proposition 5.51). 

In this paper, in section El we supply new estimates for recursively defined 
norms depending on gaugepotentials. This is perhaps our main theoretical nov¬ 
elty. The numerical method is defined and analysed in section [3] Finally inE] 
we include scalar potentials and spin, as in the Pauli equation. 
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2 Comparisons of recursive norms 

In this section we study functions defined on simplexes. Given a simplex T, we 
consider functions in H^(T) with the additional property that the restriction to 
any face T' of T is in H^(T). We denote by Hjg(.(T) the space of such functions, 
to point out the recursive nature of the construction. One advantage with this 
space with respect to hnite element analysis is that the nodal interpolator (which 
requires taking vertex values) is well defined independently of space dimension. 
This reduces the need for appeals to regularity theorems and Sobolev injections. 

Our goal is to relate different norms of such functions, in particular norms 
depending on a choice of gauge potential A, defined on T as well as on its 
faces. For simplicity we consider only the case of gauge potentials which are 
constant, so that the differential operators /S.a have constant coefficients. This 
has the advantage of guaranteeing for instance that the kernel of is a one¬ 
dimensional space of functions. 

In what follows A will denote a set consisting of constant one-forms (gauge 
potentials) attached to T and its faces, which is bounded say with respect to 
the L°° norm. In this section, most often one should think of T as a “reference” 
simplex, of diameter of order 1. If T is a “physical” simplex of diameter h and 
we map it back to a reference simplex T and pull back a gauge potential A on 
T ioT, A gets multiplied by h. Thus, if A is bounded on S', the pullbacks of A 
to reference elements will indeed live in a bounded set. Notice that the set of 
constant one-forms on each subsimplex of a given simplex is finite dimensional, 
therefore A will be compact, say in L“ norm. 

Notationwise we use T' <3 T to denote that T' is a subsimplex of T, T' <\.T 
to say that T' is a subsimplex of T which is not a point, and T' <\' T to say 
that T' is a subsimplex of T distinct from T. Given a simplicial complex T, 
denotes the set of siinplices in T of dimension k. Similarly denotes the set 
of vertices of a simplex T. 

Bounds on norms. 

Proposition 2.1. Let T be a Lipschitz domain. 

We have for u S H^(T) and A G A: 

I|w||l2(t) I|Vau||l2(t) + l|u||L2(aT)i (25) 

Proof. We consider first the case A = Q. Suppose the inequality does not hold. 
Let Un be a sequence in H^(T) such that: 

||Wn||L2(T) = 1 (26) 

||Vu„||l2(t) + ||M„||L2(aT) —t 0. (27) 

We may extract a subsequence which converges in L^(T). Since Vit„ converges 
in L^(T), this subsequence is Gauchy in H^(T), hence converges. The limit is a 
constant function which is 0 on the boundary. But it should also have norm 
equal to 1. Since this is impossible, the inequality must hold. 

Consider now the possibility that A ^ Q. Recall Kato’s pointwise a.e. in¬ 
equality: 

|V|u|(a:)| < |VAM(a:)|. 


Then one applies (051) to |m|. 
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(28) 

□ 


By induction on dimension we deduce: 

Proposition 2.2. Let T he a simplex. We have, for u £ Hjgj.(r) and A £ A: 

T'<3.T T'<I.T x£T« 

Bounds on norms. 

Proposition 2.3. Let T be a Lipschitz domain. 

We have for u £ H^(T) such that u\dT £ H^(dT) and A £ A: 

I|u||hi(t) II^AWII H-i(T) + l|u||Hi(ar)- (30) 

Proof. Suppose not. Choose sequences Un and An £ A such that: 

l|Wn||Hi(T) = 1, (31) 

I|Aa nMnlln-qT) + ll'WnllHdaT) 0- (32) 

We may suppose that An —>■ A, since A is compact. 

Let Vn be an extension to T of Un\dT, converging to 0 in H^(T). Define 
Wn = Un — Vn & Ho(’r)- We have: 


ll'a'nllnqT) 1) (33) 

IIAa„w;„||h-i(t) —t 0. (34) 

Extract a subsequence of Wn which converges weakly in H^(T) to, say, w £ 
HJ(T). We have A^rc = 0 so re = 0. So converges to 0, strongly in L^(T). 
We can also write: 

l|VA„1«n||L2(T) < l|AA„t«n||H-i(r)lkn||Hi(T) 0. (35) 

From which we deduce: 

IIVrCnIlLS ^ ||VA„ri;„||L2 + ||Wn||L2 0 

This contradicts (IHHll . 

By induction we deduce: 

Proposition 2.4. We have estimates, for u £ Hjg(.(T) and £ .4 
llwlInqT) ^ X! I|Aam||h-i(t') + |u(a;)|. 

T'<i.T T'<I.T xGTO 

Bounds on norms 

Lemma 2.5. Any element of £ H^gg(5T) has an extension to T which is in 
H^(T). More precisely, a linear bounded extension operator H^gg(5r) H^(r) 
can he constructed. 


(36) 

□ 


(37) 
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Proof. The existence of the extension can be shown exactly as in the proof of 
Proposition 3.3 in the preprint of [3]. 

The first step is to show that, for any face T' of T, if z; € H^(T') n Hj(r') 
it can be extended first to an function on the affine space spanned by T' 
and then extended to all of T by pullback and cut-off, in such a way that the 
extension is 0 on all faces of T with dimension dimT', except of course T' itself. 

The second step is a recursive construction of the extensions, detailed in 
Proposition 2.2 in m- □ 

Proposition 2.6. Let T be a simplex. For u € H^g(.(T) and A € A we have an 
estimate: 

I|w||h2(t) ||Aau||l2(t) + ll'«||H2^„(ar)- (38) 

Proof. Suppose now that the estimate does not hold. Choose a sequence Un in 
H,^gg(T) such that: 


l|Mn||H2(T) = 1, (39) 

l|Ayi„'u„||L2(T) + ||M„||H2^jaT) 0- (40) 

Using the preceding proposition, choose an extension Vn of Un\dT to T such 
that: 

lkn||H2(T) 0. (41) 

Setting Wn = Un — Vn G H^(r) n Hq(T) we have: 

l|w'n||H2(T) 1, 

We have: 


< \\AA^Wn\\h^(T) l|w^ra||L2(T) 0, (44) 

Hence Wn converges to 0 in H^(r). We then have: 

IIAzrn||L2(T) <11 Aa„'u;„||l2(t)+ (45) 

(2||24„||l= -I- ||^n||L“ + II divH„|| l°=)||u'„||hi(t), (46) 

^0. (47) 


By elliptic regularity we deduce that converges to 0 in H^(T). This 
contradicts (H3]) . 

□ 

Remark 2.1. ProDOsition l2. di does not hold on arbitrary cells. It fails for instance 
for cells that have several adjacent faces which are coplanar. 

Remark 2.2. Lemma 12.51 and its proof can be extended verbatim to for 
k > 2. However we only need it in conjunction with Proposition 12.61 which is 
limited to the case fc = 2, since elliptic regularity on simplexes does not give H^ 
estimates. 

By induction we deduce: 


(42) 

(43) 
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Proposition 2.7. Suppose u G Hjgj,(r) satisfies, for each T'<i.T, Au G L^(T'). 
Then u G H^gj,(T) and we have an estimate: 

T'a.T T'<i.T x&TO 


More bounds. In this paragraph we are interested in bounding various quan¬ 
tities by a recursive norm based only on the covariant Laplacian and norms 
attached to the whole simplex T. 

Proposition 2.8. For u G Hjgj.(T) and A € A, we have an estimate: 


y;Kx)N i: II^.4 m||h-i(T') + lkllL2(T)- 

X T'O.T 

Proof. If not, choose so that: 

\unix)\ = 1 , 

X 

X] ll^^n'“llH-i(T') + I|m||l 2(T) -t 0. 

T'-d.T 


(49) 


(50) 

(51) 


The sequence is bounded in H^gj,(r), by ProDOsition l2.4l Extract a subsequence 
which converges weakly in Hjgj.(T). The limit must be 0. This contradicts 
dsni), because, on edges, the trace operator is compact from H^(T') to vertex 
values. □ 


Proposition 2.9. Let T be a simplex. 

We have the following estimate. For all A G A and u G H):gj.(r) we have: 


X X II^^“IIh-i(T') + I|m||l2(t)- (52) 

T’CT T'C.T 

Proof. We simply write, for any subsimplex S' of T: 

I|u||l2(S) X II^AmIIh-i(T') + X 

T'<i.S xeSO 

^ X ii^^“IIh-i(t') + X 

T'd.T x^T° 

^ X II^^“IIh-i(t') + I|m||l2(t), (55) 

T'C.T 

using Propositions 12.41 and 12.81 □ 


Proposition 2.10. Let T be a simplex. For u G H):gg(T) and A G A we have a 
bound: 

X X II^^“IIh-i(T') + I|Vu||l2(t) + ||^||l°°||m||l2(t)- (56) 

T'<I.T T'a.T 
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Proof. Suppose that Un G Hjgj.(T) and An G A are sequences such that: 


^ ||Vu||l2(t') = 1, (57) 

T'd.T 

ll‘^A„'Wn||H-i(T') + + ||A„||l~ ||'«n||L2(T) 0. (58) 

T'a.T 

Suppose first that: 

liminf ||A„||l~ > 0. (59) 

n 

Then u„ converges to 0 in L^(r), and we get a contradiction using Propositions 
12.81 and [2.41 

Suppose, on the other hand, that we can extract a subsequence such that 


||24„||l°° —0. Let Un be the average of on T. 

We then have: 

llwn — 'U„||h1(T) 0. (60) 

We also have, for any face T' of T: 

||AA„(Un — '«n)||H-i(T') < II || h-i (T') + II div^ (A{t„) ||h-i(T') : (61) 

II^A^MrallH-qT') + ||24n'Un||L2(T') (62) 

II^A„M„||h- 1(T') + ll^nllu” ||firi||L2(T), (63) 

^ 0. (64) 

It follows that Un — Un converges to 0 in Hjg^(T), which contradicts (157)) . 

This concludes the proof. □ 


Proposition 2.11. Let T be a simplex. For u G H^g(,(T) and A G A, we have 
an estimate: 


|u|h2(T') 'Y2 ll^^■“llL^(T) + ll^llL“l|57u||L2(T) + ll^llL“ll'«llL2(T)• (65) 

T'<|T T'O.T 

Proof. We proceed as in the preceding proof. 

Suppose that Un G Hjgj.(T) and An G A are sequences such that: 



H2(T') — ij 

(66) 


T'<I.T 


L2(T') 

+ 1 L“ 1 Vm„| l2(T) + ll^ra L“ "“n L2(T) ^ 0. 

(67) 

T'<I.T 



Suppose first that: 

liminf > 0. 

(68) 


Then we get a contradiction from Proposition [7151 

Suppose, on the other hand, that we can extract a subsequence such that 
||^n||L°“ 0. Let Un be the affine function which best approximates in 

norm. We have: 

\\Un — Un||H2(T) ^ |Wn|H2(T) < 1- (69) 
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We also have: 


||Ayi„('u„ - £t„)||L2(T') (70) 

^ II AA„'«n||L2(T') + l|Ayl„'Un||L2(T'); (71) 

< II AA„Mn||L2(T') + ll^lnllL” II V{t„||L2(T') + l|24„ || ^<=0 || Hl^ (T') > (72) 

< II AA„Mn||L2(T') + ||24„||l“ II Vu„||l2(T) + ll^ln ||l“ llT^n II L2(T) , (73) 

^ 0. (74) 


It follows that Un — Un is bounded in H^g^(T). Extract a weakly converging 
subsequence. The limit must be harmonic on every subsimplex, hence affine. It 
is also orthogonal to affine functions, hence 0. Therefore the vertex values of 
Un — Un converge to 0. Combined with (EH) this shows that — Un converges 
to 0 in H^g(.(T), contradicting (l66)) . 

□ 


3 Numerical method 

3.1 Whitney forms 

Given a simplicial mesh, the barycentric coordinate map associated with vertex 
i is denoted A^. It is continuous, piecewise affine, has value 1 at vertex i and 0 
at the other ones, and is uniquely determined by these properties. 

We assume that all simplices in our mesh have been oriented. Given a 
simplex T of dimension k, with vertices numbered from 0 to k, the associated 
Whitney form is: 


k 

Xt = fc! ^(-l)^AjdAo A ... dAj... A dAfc. (75) 

j=o 

This sum depends on the numbering of the vertices only up to a sign, which is 
1 iff the numbering is compatible with the orientation of the simplex. 

The span of the forms Xt with T ranging through the A:-dimensional sim- 
plexes in some simplicial complex T, is denoted 211^ (T): 

2U'=(r) =span{AT : T (76) 

Such forms are also called Whitney forms. 

For instance given an edge e, oriented from the vertex e to the vertex e, the 
associated Whitney one-form is: 

Ae = AgdAe — AedAe. (77) 

When A is a general Whitney one-form, we may write: 

Al = ^AleAe, (78) 

e 

with: 

4le = [a, (79) 

J e 


10 


given that for each edge e, an orientation has been chosen. 

We will also use the following notations. If x and y are two vertices we put: 

{ Af^ if a; = e and y = e, 

—Ag ii X = e and y = e, (80) 

0 other cases. 

In fact this coincides with the definition: 

A^y= [ A, (81) 

J[xy] 

where [xy] is the oriented edge connecting x to y. 

3.2 Parallel transport 

Suppose [xy] is the oriented edge connecting x to y, and that A is a constant 
one form on it. Given a value Uy associated to y we may solve, on [xy], the 
ordinary differential equation for u: 



V.4M = 0. 

(82) 

We then have 

a relation between the vertex values of u: 



II 

(83) 

with: 

Uxy — 

(84) 


given the identity (ISTl) . 

We refer to Uxy as the parallell transport from y to x. We may notice: 


Uyx = U^. (85) 

We also use the convention: 

Uxx = l- ( 86 ) 

Discrete gauge transformations are defined as follows. We assume we have 
a value € K. associated with each vertex x G 7^. It acts on u defined at 
vertices x, with vertex values Ux, by the transformation: 

Ux exp{iax)ux- (87) 

We will always consider that gauge potentials are in 7B^{T) and that gauge 
transformations are defined by functions a G so that grad a G 2U^(T). 

However the wave functions u will be reconstructed from their vertex values, 
in several different ways, depending on a constant gauge potential A on each 
simplex T (of every dimension). These gauge potentials can in principle be 
chosen independently of each other. However, in practice, on any simplex T, 
the associated A will be the average of some globally compatible A G 211^ (T). 
For instance, we may then write, when T is a tetrahedron, for x G T: 

A{x) = A+{1/2)B X {x-xt), (88) 


II 


where B = curl A and xt is the isobarycentre of T. 

We denote by Xt [A] the space of complex valued functions it on T such that 
for every subsimplex T', we have: 


A^it = 0. (89) 

For instance ^^[O] is simply the space of affine functions. Elements of 
are uniquely determined by their vertex values, so that we may consider the 
nodal interpolator onto 

Given regular meshes Th of size h, the associated space of scalar functions 
will be denoted (note that A depends on h). 

Passing between X;,[0] and X;i[X] does not produce big errors. In fact we 
have the following estimates: 

Proposition 3.1. Let II/i denote the nodal interpolator onto X^t[0]. For u G 
X[A] we have estimates: 


||u - n,^it||L2 ^/i^||u||hi, (90) 

||Vu - Vn^,it|lL2 ^ /i||it||Hi- (91) 

Proof. Consider a simplex T of diameter 1. We write, using a Bramble-Hilbert 
type argument, followed by Proposition l2.11l 

||it — n^ii||L2(T) (92) 

T'<i.T 

4 ||A||l~||Vii||l2(t) + II^IIl“II«IIl2(t)- (93) 

Then estimate (1901) follows by scaling. 

The other estimate is proved similarly. □ 


In the other direction we note: 

Proposition 3.2. Let II/j denote the nodal interpolator onto X/i[X]. For u S 


X^i[0] we have estimates: 

||ii-nftM||L2 ^/i^lliillni, (94) 

||Vii - Vn/jit||L2 ^/i||it||Hi- (95) 

Proof. Consider a simplex T of diameter 1. We apply Proposition l2.4l 

\\u — nu||L2(T) E l|Ai-l|y|r.|. (96) 

T'<\.T 

4 E PI|l~||Vu||l.(to + P||^.||ii||l2(to, (97) 

T'<\T 

4 II^I|l°“||Vm||l2(t) + II^IIl“II'w|Il2(t)- (98) 

Then (IMll follows by scaling. 

The other estimate is proved similarly. □ 


Remark 3.1. From this last Proposition estimates of best approximation in 
can be deduced from those in X/,[0], which are well known. 
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3.3 Covariant mass matrix 

We now wish to define scalar products of functions, given their vertex values. 

For a tetrahedron T, we denote by M{T) the mass matrix on 211°(T) with 
respect to the canonical basis: 

M,y{T) = J^KXy. (99) 

This matrix is real, symmetric and positive definite. The global mass matrix, 
which may be assembled from the local mass matrices defined above, will be 
denoted M, so that: 

Mxy = J XxXy. (100) 

Definition 3.1. Given parallel transports U between neighbouring vertices 
(that is, those connected by an edge), subject to (1^ and (IMl) . we define a 
covariant -product, as follows. Given complex scalar fields u and v, with well 
defined vertex values, we set: 

{u, u){7 — ^ ^ '^Ic^xyUxyVy, (101) 

xy 

where the matrix M is the mass matrix already defined in (IIOOII . 


This matrix may also be assembled from terms local to each tetrahedron. 


Proposition 3.3. This covariant scalar product S101\) 
invariant, under transformations : 

is hermitian and gauge 

Ux exp{iax)ux, 

(102) 

Uxy exp{iax)Uxy exp{-iay), 

(103) 

Vy 1—>■ exp(iay)uy. 

(104) 

Proof, (i) Hermitian: 


(l?, u)lJ — ^ ^ '^x^xyUxyUy-) 

(105) 

xy 


II 

M 

(106) 

xy 


= ('U,v)u- 

(107) 

(ii) Gauge invariance is trivial. 

□ 


Proposition 3.4. The covariant scalar product fMD is h'^-conforming with 
respect to the norm on the space X;i[0] associated with a regular mesh of 
width h. In other words, for u,v G we have: 


\{u,v) - {u,v)u\ 4 ^■^||M||Hi|kllHi- (108) 

Proof. We work on a simplex T of diameter h. 

We want to estimate the error: 

{u,v) - {u,v)u, (109) 
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And we do this by decomposing u and v according to: 

u(x) = ut + u'rp • (x — xt), with Uj- = gradu, (HO) 

v{x) = vt + v't ' (x — Xt), with u't = grad?;, (HI) 

where ut denotes the value of u at the isobarycenter of T, which is denoted xt, 
and u'rp is the gradient of u. Likewise for v. 

We now treat the four terms this decomposition gives: 

(i) We have: 


{ut,vt) — {ut,vt)u = — 1)ut, (H2) 

xy 

= ulpMa,y (C0S{ Aa,y) -1)VT, (H3) 

xy 

4 ft,^||u||L2(T)||w||L2(T)- (H4) 

The trick was to symmetrize in x and y, use (18511 and | cos(A 2 ,y) — 1| ^ 

(ii) We have: 

{uT,v'j, ■ (jj - yr)) - {uT,v'rp ■ {y - yT))u, (H5) 

= ^u)pM,,y{U^y -1){vt ■ {y-yr)), (H6) 

xy 

^/i^||u||l2(t)|| grad?;||L2(T)- (H7) 

Here we combined \Uxy — 1| ^ /i and \y — yrl < h. 

(iii) There is a similar term where u and v exchange roles. 

(iv) Finally, the last term, involving two gradients, yields a factor h^. 

This completes the proof. □ 


Proposition 3.5. The covariant scalar product ([M]) is -conforming with 
respect to the norm on the space A/i[A] associated with a regular mesh of 


width h. In other words, for u,v G 

|(m,v) - {u,v)u\ =4 /i^||M||Hi|kllHi- (H8) 

Proof. We let H^ denote the nodal interpolator onto Xh [0]. We have, by Propo¬ 
sition 13.11 

\{u,v) - {IlhU,Ilhv)\ 4 ^^IImIIhiIIwIIhi, (H9) 

and, by Proposition [SHI 

|(n,iU,n;,?;) - {T].hU,I\hv)u\ ^ ft.^||n;iM||Hi||n,j?;||Hi. (120) 

We conclude by the stability of H/j : Xh [A] —>■ Xh [0] in norm, which can be 
deduced from Proposition 13.II □ 

3.4 Covariant stiffness matrix 

Lemma 3.6. La A be a constant one-form on the edge [xy\. We consider 
functions u : [xy] —>■ C satisfying Aau = 0. Then: 

{VAu{x)){y -x) = u{y) exp,{iAxy) - u{x), (121) 

{VAu{y)){y -x) = u{y) - exp{-iAxy)u{x). (122) 
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Proof. We parametrize the edge linearly from 0 to 1, with a variable t. The 


solutions to Xau 

= 0 have the form: 



u{f) = (a -I- bt) exp(—iAajyt). 

(123) 

We get: 

u(0) = a. 

(124) 


u(l) = {a + h) eyzpi-iAxy). 

(125) 

We also compute: 

(V -\- iA)u{t) = eyzp{—iAxyt)b. 

(126) 

We deduce: 

(V-kiA)u(O) =6, 

(127) 


= u(l) eyiY>{iAxy) - u(0). 

(128) 

and likewise: 

(V -1- iA)u{\) = eyzp{—iAxy)b, 

(129) 


= u(l) - exp(-iA2,y)M(0). 

(130) 

From this the lemma follows. 

□ 


We use these identities as follows. Given a mesh and a Whitney one-form A 
on it, we define, on every simplex. A, to be the average of A. We remark that 
on edges A = A. Associated with A we have the space X[A] of scalar functions, 
on the mesh, determined by imposing = 0 on every simplex (of every 

dimension). Then for u £ ^[A], the vertex values of Xau are computable, using 
(1121111221) . Thus we may interpolate it onto affine vector fields in a computable 
way. We now detail how this leads to a numerical method. 

Let r be a tetrahedron. We consider, at each vertex x and for each edge 
[xy] emanating from x, the tangent vector Txy = y — x. At any vertex a;, the 
three tangent vectors, pointing to the three other vertices of the tetrahedron, 
constitute a basis of K^. We denote by y,xy the dual basis, i.e.: 


l^xzij'xy') — 


1 \iy = z, 
0 ii y ^ z. 


(131) 


Given nodal values for a scalar function u on T, we construct an affine one- 
form V by setting for each pair of vertices (x, y): 


Vxy = v{x){Txy) = UxyU{y) - u{x) . 

and, summing over pairs of vertices we define the vectorfield: 

^ ^ ^ '^xy^xl^xy- 


(132) 


(133) 


By Lemma 1531 if u £ A[A], u is then the affine vectorfield on T coinciding with 
Vau at the vertices. 


15 




The scalar product on affine one forms (or vector fields) can be expressed 
with the scalar mass matrix, as follows. Let v and v' be affine one forms. Define 
the numbers v^y by (11321) . and proceed similarly for and v'. 

v'^ -v' = ^ vlyv'^t (134) 

xy,zt 

The covariant scalar product on affine one forms is defined by setting; 

{v,V )lJ = ^ ( '^xyUxz'Vzt ■ fJ-ztMxz- (135) 

xy,zt 

Definition 3.2. Let u and v be functions with well defined vertex values. We 
define a modified bilinear form, as a sum over tetrahedra: 

a{u,v) = E driu^v)^ (136) 

T 

where the contribution of tetrahedron T is: 

^ ^ {Uxy'^y ’^x)^Uxz{Uzt'^t '^z^f^xy ' f^zt^xz{T^. (137) 

xy,zt 

A crucial identity, which follows from Lemma iTbl is that for u,v € Ai[A]: 

a(M, w) = (IIVau, LIVau);/. (138) 

The interpolation operator 11 appearing here, is the nodal interpolation onto 
affine one forms. 

Proposition 3.7. The discrete stiffness matrix defined by JlgT] ) is gauge in¬ 
variant (in the same sense as in Proposition roi) and hermitian. 

Proof. Trivial. □ 

Next we examine the consistancy of the discrete stiffness matrix. The key 
estimate is the following one: 

Proposition 3.8. Let II/i denote interpolation onto affine one forms. Choose 
u € Xh[A\. We have : 

||Vau - n;iVAu||L2 h^||u||Hi. (139) 

Proof. We work first on a simplex T of diameter one. We have an estimate: 

||Vau — n/jVA'«||L2(T) ^ ^ ||AVau||l2(t')- (140) 

T'O.T 

Consider now a face T' of T. We write: 

AVAU = • VVau. (141) 

We consider the three terms on the right hand side. 
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(i) First term: 


+ A^(iylu), (142) 

= A^(*Aw), (143) 

since u € A[j 4]. Suppose more generally that v is any affine function, and that 
u is arbitrary. Then: 

A^{vu) = v{A^u) + 2Vv ■ Vu + 2i{A ■ Vv)u. (144) 

Letting the role of v be played by A, we may continue from (1143^ : 

||A^(*Au)||l2 ^ ||B||l.»||Vu||l 2 + ||i||L~||i3 ||L~||u||l=. (145) 

(ii) Second term: 

II I^I^VauIIl. < ||i||2„„||Vu||L2 + ||i||^..||A||L~||u||L2. (146) 

(hi) Third term: 


I|i-VVAU||L2 ^ ||i||L~||heSSu||L2 + ||i||L~||S||L-||w||L^ + P||L~||4l||L=||Vu||L2. 

(147) 

Summing the three terms, we get, on T' : 


||AVau||l 2 <||dl||L~||hessu||L2 + (P||2„„ + ||B||Lo=)||Vu||L2+ (148) 

(||A||3^ + P||L»||i3||L»)||u||L2. (149) 

Using Propositions 12.9112.10112.111 the right hand side terms on T' may be 
bounded by terms attached to T, as follows: 

IIVau - nVAu||L2(T) 4{\\A\\l^ + ||il||L<»)||Vu||L2(T)+ (150) 

(PII^o= + I|24||l==||b||l==)IHIl^(t). (151) 


Finally one concludes by scaling, noting that A scales like a one-form and B 
like a two-form. 

□ 


We also notice the variant: 

Proposition 3.9. Let Il/i denote interpolation onto affine one forms. Choose 
u £ Xh [A]. We have : 


||VVau - VU/iVamIIl^ < li||M||Hi- (152) 

Proof. We go through the preceding proof to obtain the following variant of 

(fTCUD : 

II VVau - Vn;,VAu||L2(r) ^(Pll^o. + ||B||l-)||Vu||l2(t)+ (153) 

(PllE== + l|24||L=o||i?||L~)||«||L=(T). (154) 

Then the scaling gives the factor h this time. □ 
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Proposition 3.10. For u,v in X/i[j4] we have an estimate: 

I J{yAu)^-yAV - a(u,v)\ ^ h^\\u\\-iii\\v\\-iii. (155) 

Proof. We first use Proposition 13.81 to get: 

I y'(VAu)t.V^u-{n,,VAU,nVAu)H/i2||u||Hi||u||Hi. (156) 

Then we use Proposition [3il] (generalised from scalar to vector fields): 

|(n/,VAU,n,iVAu) - (n^VAU,n?,VAu)(7| ^ /i^||n?iVAu||Hi||n^VAu||Hi- (i57) 
The rest of the proof is devoted to bounding the right hand side of this estimate. 


Using Proposition 13.81 we get: 

||n?iV^u||L2(T) ||V^m||l2(t) + ^^||u||hi(t), (158) 

^ I|u||hi(t)- (159) 

Using Proposition 13.91 we get: 

||Vn^V^'u||L2(T) IIVV^m||l2(t) + 1i||^^IIhi(t)- (160) 

On a simplex of diameter 1 we have, using Proposition 12.Ill 

||VVau||l2(t) I|VVu||l2(t) + II^IIIl”IIVu||l2(t) + I|11||l“||u||l2(t); (161) 

^ P||l~||Vu||l2(t) + (Pll^. + ||i3||L~)||u||L2(T). (162) 

Since the left and the right hand side scale the same way, we deduce: 

IIVn/iVA'u||L2(T) II'w||hi(t)- (163) 

From (|159|1 and (116311 we may conclude: 

||II;jVau||h1(T) I|w||h1(T)- (164) 

Inserting this in (11571) completes the proof. □ 

3.5 Conclusions 

Summing up, the situation is as follows. 


Consider the eigenvalue problem: Find u G Hj(S') and F G R such that for 
alluGHi(S'): 

a(u,v) = F(u,v}. (165) 

Since we assume that the domain S is convex and that the gauge potential A is 
smooth, elliptic regularity holds, in the sense that the solution operator for 
maps L2(S') to Hi(S') C 

We do the Galerkin formulation on the space A/j[A] attached to a mesh of 
width h. The order of convergence for discrete eigenvectors is h in norm and 
hf in lA norm, since these are the orders of best approximation, see Remark l3. II 
The order of convergence of the eigenvalue is deduced to be (see in particular 
Lemma 3.1 in CD- 
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Then we consider the modified formulation: Find u G and E such 

that for all v G Xh [A]: 

d{u,v) = E{u,v)u, ( 166 ) 

where the modified bilinear form d was defined in Definition 13.21 whereas (•, ■)u 
was defined in Definition 13.11 These modifications produce an error of order 
in norm, as was shown in ProDOsitions l3.10l and 13.51 Therefore the preceding 
orders of convergence for eigenvectors and eigenvalues are maintained (for the 
eigenvalue see Lemma 5.1 in 0)- 

Finally we may interpolate the eigenvector u G onto X/i[0] and still 

get the same orders of convergence, using Lemma [3.21 

Of course we could also consider that we do the variational formulation of 
(11651) on (rather than and that the discretization (11661) . where d is 

defined explicitely from vertex values in (11371) . constitutes a variational crime on 
However for elements u,v € ^^,[0] the formula (11381) . which is essential 
to our analysis, is then no longer true. In fact we expect the error of consistency 
in norm to be h in this interpretation, whereas it was h? in the preceding 
one. This would ruin the analysis, even though we have described the same 
numerical method. 

Finally it must be noted that the above analysis requires H to be a Whitney 
form on each grid. There is an implicit step where, given H on S' one first 
approximates it on the grid Th, to be able to use the previous analysis. In 
general this step produces an error of order h for smooth A, jeopardizing the 
above analysis. For details on how order h can be obtained, under weaker 
hypotheses, see mm- But in the important case of a constant magnetic field 
B, the magnetic vectorpotential has (globally) the form: 

A{x)=Ao + il/2)Bxx, (167) 

for some constant Aq . Then the approximation of A by Whitney forms is exact 
and our analysis applies. We consider this case to be sufficiently important to 
justify our new method. 


4 Extension to the Pauli equation 


A Pauli wave function is a complex valued two-component spinor tjj in Hj(S') 0 
C^. We write: 


V' = 


'00 

'01 


(168) 


Let a = (cTi,(72,(Ta) be the hermitian and unitary Pauli matrices collected in a 
vector. The components are: 


(Tl = 


0 

1 




t73 = 


1 

0 



In natural units the time-dependent Pauli equation reads: 


- ^(i?-V a)^0 


idvtpj 


(169) 


(170) 
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where S/a'^P = grad^' + iAip is the covariant spatial gradient of and dy'p = 
ip + iVip is the covariant time derivative. If we assume that the time-dependence 
is tpixjt) = where E gR, then we get the Pauli eigenvalue equation 


-^{a-VA)^'ip + Vi; = EP;. 


(171) 


From the identity 


-f 'i ^ 


(172) 


where e is the totally antisymmetric Levi-Civita symbol with £123 = 1, we can 
rewrite the equation as 


— + S ■ cml A)^!} + VIp = Eip. (173) 

The variational formulation of the Pauli eigenvalue problem consists in find¬ 
ing Ip G {E[q{S) 0 C^) and E G R, tp ^ 0, such that for all (p G H),(S') 0 

a{ip, (p) + bpip, (p) + c{ip, (p) = E{ip, (p), (174) 

where a(-, •), &(•, •), and c(-, •) are the bilinear forms given by: 

a('>P,P) = (VaiP,VaP}, 

= (V'P,P), (175) 

c{ip,p) = -{[S ■ cnT\A)ip,(p). 

Equation (11741) remains invariant under gauge transformations (ITSl [HI) . We 
proceed to define gauge-invariant discretizations. 

The bilinear forms a and (•,•) are discretized as before, and the analysis 
carries over straighforwardly to spinors. 

The forms b and c are treated in analogy with the previously defined covari¬ 
ant scalar product. 

Explicitely we define: 

H'fp, (P) = ■ UxyPy (176) 

xy -I 


If E is a smooth potential on S one can approximate it by piecewise linears, to 
second order, before plugging it in to (11761) . as was done in [7]. 

We also define: 



Here we rely on the fact that a • curlH is a constant matrix on each T. 

Since there are no novelties in the proofs, we just state the conclusions: 
The above discretization technique yields second order convergence under the 
same hypotheses as before, which include constant magnetic fields and smooth 
scalar potentials. Lower order convergence rates can be obtained under weaker 
hypotheseses, as in [7]. 
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